/**************************************************************************//**
 * @file    min_search_parabolic_interpolation.c
 * @author  <a href="http://mymathlib.webtrellis.net">Mathematics Source Library</a>
 * @date    02-08-2010
 *
 * <b>Description</b>
 *
 * The extremum of the parabola \f$ y = cx2 + bx + a \f$ where \f$ c != 0 \f$ is located at \f$ x* = - b / 2c \f$ . If
 * \f$ c < 0 \f$ , then the extremum is a maximum and if \f$ c > 0 \f$ , then the extremum is a minimum.\n
 *
 * Suppose that we are given a smooth function \f$ f : [α,β] → R \f$ , where \f$ α < β \f$ for which we seek to
 * locate a relative minimum. Assume that also we are given the initial data \f$ (x1, f (x1) ) \f$ ,
 * \f$ (x2, f (x2) ) \f$ , and \f$ (x3, f (x3) ) \f$ , where \f$ α ≤ x1 < x2 < x3 ≤ β \f$ such that \f$ f (x1) > f (x2) < f( x3) \f$ .
 * Then the location of minimum, x*, of the parabola through the same three points must lie between
 * x1 and x3. If \f$ x* = x2 \f$ then perturb it by a small amount \f$ η \f$ either towards x1 or x3 depending on
 * whether \f$ f (x1) \f$ is less than or \f$ f (x3) \f$ or not. Then either \f$ x* < x2 \f$ or \f$ x* > x2 \f$ . \f$ If x* < x2 \f$ then
 * either \f$ f (x*) ≤ f (x2) \f$ or \f$ f (x*) > f (x2) \f$ . If \f$ f (x*) ≤ f (x2) \f$ then set \f$ x3 ← x2, x2 ← x* \f$ or if
 * \f$ f (x*) > f (x2) \f$ then set \f$ x1 ← x* \f$ . The process is analogous to that above if \f$ x* > x2 \f$ . Repeat the
 * process unless the length of the interval of uncertainty is less than a preassigned tolerance in
 * which case terminate the procedure.
 *
 * <b>Function List</b>
 *
 * int \ref Min_Search_Parabolic_Interpolation( double (*f)(double), double* x0, double* fx0, double tolerance )
 *
 * Given initial data ( x0[0], fx0[0] ), ( x0[1], fx0[1] ), and ( x0[2], fx0[2] ) where
 * x0[0] < x0[1] < x0[2], fx0[0] > fx0[1] < fx0[2] and fx0[i] = f ( x0[i] ), i = 0, 1, 2
 * this routine reduces the interval of uncertainty [ x0[0], x0[2] ] by an extension of the
 * iterative process described above. The iteration halts when the length of the interval of
 * uncertainty relative to its midpoint is less than tolerance. Upon return, x0[0] is the lower
 * endpoint of the interval of uncertainty, x0[2] is the upper endpoint and x0[1] is the current
 * best estimate of the location of the minimum. The array fx0 will contain the value of f ( x )
 * evaluated at the three points.
 *
 * int <b>Max_Search_Parabolic_Interpolation</b>( double (*f)(double), double* x0, double* fx0, double tolerance )
 *
 * Given initial data ( x0[0], fx0[0] ), ( x0[1], fx0[1] ), and ( x0[2], fx0[2] ) where
 * x0[0] < x0[1] < x0[2], fx0[0] < fx0[1] > fx0[2] and fx0[i] = f ( x0[i] ), i = 0, 1, 2
 * this routine reduces the interval of uncertainty [ x0[0], x0[2] ] by an extension of the
 * iterative process described above. The iteration halts when the length of the interval of
 * uncertainty relative to its midpoint is less than tolerance. Upon return, x0[0] is the lower
 * endpoint of the interval of uncertainty, x0[2] is the upper endpoint and x0[1] is the current
 * best estimate of the location of the maximum. The array fx0 will contain the value of f ( x )
 * evaluated at the three points.
 *
 * @see http://mymathlib.webtrellis.net/c_source/optimization/nonlinear/one_dim/min_search_parabolic_interpolation.c
 * @see http://mymathlib.webtrellis.net/optimization/nonlinear/one_dim/parabolic_interpolation.html
 *
 *****************************************************************************/

#include <float.h>                          // required for DBL_EPSILON
#include <math.h>                           // required for fabs(), sqrt()

#define SQRT5 2.236067977499789696

static const double alpha = 0.5 * (3.0 - SQRT5);

//                        Internally Defined Routines

static int Stopping_Rule(double x[], double tolerance);
static int Parabolic_Minimum( double x[], double y[], double *xmin);
static void Check_New_Test_Point( double x[], double *xmin,
                                     double min_displacement, double epsilon );
static void Left_Interior_Point( double (*fn)(double), double x[], double y[] );
static void Right_Interior_Point( double (*fn)(double), double x[], double y[] );


/*! @brief Minimum Search function by using parabolic interpolation.
  *
  *    This routine can be used to locate the minimum of a unimodal function
  *    within a preassigned tolerance or a local minimum of a multimodal
  *    function within a preassigned tolerance.  The function should not be
  *    linear on any closed and bounded subinterval with nonempty interior.
  *
  *    Initially 5 points are defined \f$ (x[i], fx[i]), i = 0,1,2,3,4 \f$ by setting
  *    \f[ x[0] = x[1] = x0[0] \f]
  *    \f[ x[2] = x0[1] \f]
  *    \f[ x[3] = x[4] = x0[2] \f] and
  *    \f[ fx[i] = f(x[i]). \f]
  *    The procedure then fits a parabola through the three interior points.
  *
  *    Given points \f$ (x1, f(x1)) \f$ , \f$ (x2, f(x2)) \f$ , \f$ (x3, f(x3)) \f$ with \f$ x1 < x2 < x3 \f$
  *    and \f$ f(x1) > f(x2) < f(x3) \f$ , the parabolic interpolation method
  *    estimates a local minimum of \f$ f(x) \f$ for \f$ x1 < x < x3 \f$ by repeated
  *    estimation of the minimum of the parabola through three points. The
  *    three points being two points which bound the minimum and the current
  *    location of the observed minimum.  The iteration is terminated when
  *    the relative length of the interval bounding the minimum is less than
  *    or equal to a preassigned tolerance.
  *
  *  <b>Example:</b>
  * @code
  *     extern double f(double);
  *     double x[] = { -1.0, 0.0, 1.0 };
  *     double fx = { f(x[0]), f(x[1]), f(x[2]) };
  *     double tolerance = 1.0e-2;
  *     int err;
  *
  *             (verify fx[1] < fx[2] and fx[1] < fx[0] otherwise
  *                                 choose different initial points)
  *
  *     err = Min_Search_Parabolic_Interpolation(f, x, fx, tolerance )
  *     if (err < 0 ) printf("Error in input data\n");
  *     else printf("local minimum is between (%12.6e,%12.6e) and
  *                               (%12.6e,%12.6e)\n", x[0],fx[0], x[2],fx[2]);
  * @endcode
  *
  * @param[in] (*f)      Pointer to a user-defined function of a real variable (type double) returning a real number (type double).
  * @param[in] *x0       On input, \f$ x0[0] \f$ is the initial lower endpoint and \f$ x0[2] \f$ is the upper
  *                      endpoint of an interval which contains a local minimum. \f$ x0[1] \f$ is
  *                      an interior point between \f$ x0[0] \f$ and \f$ x0[2] \f$ such that
  *                      \f$ f(x0[1]) < f(x0[0]) \f$ and \f$ f(x0[1]) < f(x0[2]) \f$ .  On output, \f$ x0[0] \f$ is
  *                      the lower endpoint and \f$ x0[2] \f$ is the upper endpoint of the final
  *                      subinterval which contains the local minimum. \f$ x0[1] \f$ is the current
  *                      best estimate of the location of the local minimum.
  * @param[in,out] *fx0  On input and on output, \f$ fx0[i] \f$ is the value of \f$ f(x) \f$ evaluated at \f$ x0[i], i=1,2,3 \f$ .
  * @param[in] tolerance A parameter which controls the termination of the procedure.
  *                      When the magnitude of \f$ x[3] - x[1] \f$ is less than "tolerance" *
  *                      \f$ ((x[3] + x[1]) / 2 + eta) \f$ , the procedure is terminated, where eta
  *                      is <b>1</b> if \f$ *x3 + *x1 <= 1 \f$ and <b>0</b> otherwise. <b> "tolerance' should be
  *                      is strictly positive number </b>.  If a nonpositive number is passed,
  *                      then tolerance is set to sqrt(DBL_EPSILON).
  *
  * @return Returns an integer value about succeess or fail.
  * @retval  0   Success
  * @retval -1   The three points are collinear.
  * @retval -2   The extremum of the parabola through the three points is a maximum.
  * @retval -3   The initial points failed to satisfy the condition that
  *              \f$ x1 < x2 < x3 \f$ and \f$ fx1 > fx2 < fx3 \f$ . */
int Min_Search_Parabolic_Interpolation(double (*fn)(double), double* x0,
                                               double* fx0, double tolerance ) {
	double x[5] = { x0[0], x0[0], x0[1], x0[2], x0[2] };
	double y[5] = { fx0[0], fx0[0], fx0[1], fx0[2], fx0[2] };
	//double min_displacement = 0.01 * tolerance;
	//double epsilon = 0.333333 * tolerance;
	double xmin, ymin;
	int err = 0;

		// Verify that the input data satisfy the conditions:
		// x1 < x2 < x3 and fx1 >= fx2 =< fx3 but not fx1 = fx2 = fx3.

	if ( (x[0] > x[2]) || (x[4] < x[2]) ) return -3;
	if ( (y[0] < y[2]) || (y[2] > y[4]) ) return -3;
	if ( (y[0] == y[2]) && (y[2] == y[4]) ) return -3;

		// If the input tolerance is nonpositive, set it to the default.

	if (tolerance <= 0.0) tolerance = sqrt(DBL_EPSILON);
	double min_displacement = 0.01 * tolerance;
	double epsilon = 0.333333 * tolerance;

		// Find the three interior points.

	Left_Interior_Point( fn, &x[0], &y[0] );
	if ( Stopping_Rule(x, tolerance) ) return 0;
	Right_Interior_Point( fn, x, y );

		// Iterate until the location of the minimum
		// is found within the specified tolerance.


	while (!Stopping_Rule(x, tolerance ) ) {
		err = Parabolic_Minimum( x, y, &xmin);
		if (err) break;
		Check_New_Test_Point( x, &xmin, min_displacement, epsilon);
		ymin = fn(xmin);
		if ( xmin < x[2] ) {
			if ( ymin >= y[1] ) {
				//printf("Left Interior Point Before\n");
				x[0] = xmin; y[0] = ymin;
				Left_Interior_Point( fn, &x[0], &y[0] );
				//printf("Left Interior Point\n");
			}else if ( ymin >= y[2] ) {
				x[0] = x[1]; y[0] = y[1]; x[1] = xmin; y[1] = ymin;
			}else{
				x[4] = x[3]; y[4] = y[3]; x[3] = x[2]; y[3] = y[2];
				x[2] = xmin; y[2] = ymin;
			}
		}else{
			if ( ymin >= y[3] ) {
				//printf("Right Interior Point Before\n");
				x[4] = xmin; y[4] = ymin;
				Right_Interior_Point( fn, x, y );
				//printf("Right Interior Point\n");
			}else if ( ymin >= y[2] ) {
				x[4] = x[3]; y[4] = y[3]; x[3] = xmin; y[3] = ymin;
			}else{
				x[0] = x[1]; y[0] = y[1]; x[1] = x[2]; y[1] = y[2];
				x[2] = xmin; y[2] = ymin;
			}
		}
	}

		// If no errs were detected, return the final
		// interval and an internal point together with
		// the value of the function evaluated there. durduruldu      ./siklon

	if (err == 0) {
		x0[0] = x[1];
		x0[1] = x[2];
		x0[2] = x[3];
		fx0[0] = y[1];
		fx0[1] = y[2];
		fx0[2] = y[3];
	}
	return err;
}

/*! @brief Decides stop or continue.
  *
  *     Given an ordered array \f$ x[0] < x[1] < ... < x[4] \f$ , this routine returns
  *     TRUE (1) if \f$ x[3] - x[1] < epsilon * (|xm| + eta) \f$ or FALSE (0)
  *     otherwise, where \f$ f(x[2]) \f$ is the minimum among \f$ f(x[i]), i = 0,...,4, \f$
  *     epsilon = tolerance, and eta = 1.0 if \f$ xm <= 1.0 \f$ otherwise
  *     eta = 0.0.
  *@note For C compilers which support the 'inline' type qualifier, this routine
  *      could be declared with the 'inline' qualifier.
  * @param[in] x[] array to be check.
  * @param[in] tolerance Tolerance value.
  * @retval 0  Bound is not within the specified tolerance.
  * @retval 1  Bound is within the specified tolerance. */
static int Stopping_Rule(double x[], double tolerance){
	double xmid = 0.5 * fabs( x[3] + x[1] );
	double eta = (xmid <= 1.0) ? 1.0 : fabs(xmid);
	//printf("x[3]-x[1] = %g\t tolerance*eta = %g\n",(x[3]-x[1]),(tolerance*eta));
	return ( ( x[3] - x[1] ) <= tolerance * eta ) ? 1 : 0;
}

/*! @brief Finds the parabolic minimum
  *
  *     Given three noncollinear points \f$ (x[1],y[1]), (x[2],y[2]) \f$ , and
  *     \f$ (x[3],y[3]) \f$ such that \f$ x[1] < x[2] < x[3] \f$ and \f$ y[1] > y[2] < y[3] \f$ find
  *     the location of the minimum, \a xmin, of the parabola through these
  *     points.
  *@note For C compilers which support the 'inline' type qualifier, this routine
  *      could be declared with the 'inline' qualifier.
  * @param[in] x[] array to be check.
  * @param[in] tolerance Tolerance value.
  * @retval  0  Success
  * @retval -1  The three points are collinear.
  * @retval -2  The extremum of the parabola through the three points is a maximum. */
static int Parabolic_Minimum( double x[], double y[], double *xmin) {
	double d1 = (x[3] - x[2]) * (y[1] - y[2]);
	double d2 = (x[2] - x[1]) * (y[3] - y[2]);
	double denominator = d1 + d2;
	double numerator = (x[3] + x[2]) * d1 + (x[2] + x[1]) * d2;

	if (denominator == 0.0) {                   // the points are collinear.
		*xmin = (y[1] < y[3]) ?  x[1] : x[3];
		//printf("double denominator = d1 + d2;\tdouble %.6f = %.6f + %.6f;\n", denominator, d1, d2);
		//printf("(y[1] < y[3]) ?  x[1] : x[3]\t(%.6f < %.6f) ?  %.6f : %.6f\n",y[1], y[3], x[1], x[3]);
		return -1;
	}
	if (denominator < 0.0) {                    // local maximum
		*xmin = (y[1] < y[3]) ?  x[1] : x[3];
		return -2;
	}
	*xmin = 0.5 * numerator / denominator;
	return 0;
}

/*! @brief Checks for new test point.
  *
  * This routine alters the test point, \a *xmin, in the event that \a *xmin is
  * within \a min_displacement of one of the nodes \f$ x[1] \f$ , \f$ x[2] \f$ , and \f$ x[3] \f$ .
  *
  *@note For C compilers which support the 'inline' type qualifier, this routine
  *      could be declared with the 'inline' qualifier.
  *
  * @param[in] x[] array to be check.
  * @param[out] *xmin minimum x point value
  * @param[in] min_displacement Minimum Displacement Value
  * @param[in] epsilon Tolerance value.
  * @return void */
static void Check_New_Test_Point( double x[], double *xmin,
                                     double min_displacement, double epsilon ) {
	double eta;
	double d1, d2;

				// Compare the location of the minimum of the
				// parabolic fit to the center node.

	eta = ( x[2] == 0.0 ) ? 1.0 : fabs(x[2]);
	if ( fabs(*xmin - x[2]) < (min_displacement * eta) ) {
		if ( *xmin == x[2] )
			if ( (x[3] - x[2]) > ( x[2] - x[1] ) )
				*xmin = x[2] + alpha * (x[3] - x[2]);
			else *xmin = x[2] - alpha * (x[2] - x[1]);
		else if ( *xmin < x[2] )
			if ( x[2] - x[1] < epsilon * eta )
				*xmin = x[2] + (x[2] - x[1]);
		else {
			d1 = epsilon * eta;
			d2 = alpha * (x[2] - x[1]);
			*xmin = x[2] - ( (d1 < d2) ? d1 : d2 );
		}
	else
		if ( x[3] - x[2] < epsilon * eta )
			*xmin = x[2] - (x[3] - x[2]);
		else {
			d1 = epsilon * eta;
			d2 = alpha * (x[3] - x[2]);
			*xmin = x[2] + ( (d1 < d2) ? d1 : d2 );
		}
	return;
   }

				// Compare the location of the minimum of the
				// parabolic fit to the left-most node.

	eta = ( x[1] == 0.0 ) ? 1.0 : fabs(x[1]);
	if ( (*xmin - x[1]) < (min_displacement * eta)) {
		d1 = epsilon * eta;
		d2 = alpha * (x[2] - x[1]);
		*xmin = x[1] + ( (d1 < d2) ? d1 : d2 );
		return;
	}

				// Compare the location of the minimum of the
                // parabolic fit to the right-most node.

   eta = ( x[3] == 0.0 ) ? 1.0 : fabs(x[3]);
	if ( (x[3] - *xmin) < (min_displacement * eta)) {
		d1 = epsilon * eta;
		d2 = alpha * (x[3] - x[2]);
		*xmin = x[3] - ( (d1 < d2) ? d1 : d2 );
		return;
	}
   return;
}

/*! @brief Selects new test point for left
  *
  *     Select the new test point using the points \f$ (x[0],y[0]), (x[2],y[2])  \f$ ,
  *     and \f$ (x[4],y[4]) \f$ ,  where \f$ x[0] < x[2] < x[4] \f$ and \f$ y[0] > y[2] < y[4] \f$ .
  *     Upon return the points \f$ (x[0],y[0]) \f$ , \f$ x[1],y[1]) \f$ , \f$ (x[2],y[2]) \f$ ,
  *     \f$ (x[4],y[4]) \f$ are defined with \f$ x[0] < x[1] < x[2] < x[4] \f$ and
  *     \f$ y[0] > y[1] > y[2] < y[4] \f$ , where \f$ y[i] = f(x[i]) \f$ .
  *
  * @param[in] *f control function pointer
  * @param[in] x[] x array
  * @param[in] y[] y array
  * @return void */
static void Left_Interior_Point( double (*fn)(double), double x[], double y[] ) {
	int not_found = 1;
	while (not_found) {
		//char ch = getchar();
		x[1] = x[2] - alpha * (x[2] - x[0]);
		y[1] = fn(x[1]);
		//printf("y[1]=%12.50e\ty[2]=%12.50e\n",y[1],y[2]);
		if ( y[1] >= y[0] ) {
			x[0] = x[1]; y[0] = y[1];
		if( y[1] == y[0] ){
			if(y[1] == y[2]){
				not_found = 0;
			}
		}
		//printf("y[1]=%12.50e\ty[0]=%12.50e\n",y[1],y[0]);
      	//printf("y[1] >= y[0]\n");
		}else if ( y[1] >= y[2] ) {
			not_found = 0;
			//printf("y[1] >= y[2]\n");
		}else{
			x[4] = x[3]; y[4] = y[3];
			x[3] = x[2]; y[3] = y[2];
			x[2] = x[1]; y[2] = y[1];
			//printf("else\n");
		}
   }
}

/*! @brief Selects new test point for right
  *
  *     Given four points \f$ (x[0],y[0]) \f$ , \f$ (x[1],y[1]) \f$ , \f$ (x[2],y[2]) \f$ , and
  *     \f$ (x[4],y[4]) \f$ , such that \f$ x[0] < x[1] < x[2] < x[4] \f$ and \f$ y[0] > y[1] >
  *     y[2] < y[4] \f$ find points \f$ (x[0],y[0]), (x[1],y[1]), (x[2],y[2]) \f$ ,
  *     \f$ (x[3],y[3]) \f$ and \f$ (x[4],y[4]) \f$ , such that \f[ x[0] < x[1] < x[2] < x[3] < x[4] \f]
  *     \f[ y[0] > y[1] > y[2] < y[3] < y[4] \f] , where \f$ y[i] = f(x[i]) \f$ .
  *
  * @param[in] *f control function pointer
  * @param[in] x[] x array
  * @param[in] y[] y array
  * @return void */
static void Right_Interior_Point( double (*fn)(double), double x[], double y[] ) {
	int not_found = 1;
	while (not_found) {
		//char ch = getchar();
		x[3] = x[2] + alpha * (x[4] - x[2]);
		y[3] = fn(x[3]);
		//printf("y[3]=%12.50e\ty[2]=%12.50e\n",y[3],y[2]);
		if ( y[3] >= y[4] ) {
			x[4] = x[3]; y[4] = y[3];
			if( y[3] == y[4] ){
				if(y[3] == y[2]){
					not_found = 0;
				}
			}
			//printf("y[3]=%12.50e\ty[4]=%12.50e\n",y[3],y[4]);
			//printf("y[3] >= y[4]\n");
		}else if ( y[3] >= y[2] ){
			not_found = 0;
			//printf("y[1] >= y[2]\n");
		}else {
			x[0] = x[1]; y[0] = y[1];
			x[1] = x[2]; y[1] = y[2];
			x[2] = x[3]; y[2] = y[3];
			//printf("else\n");
		}
	}
}
